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Abstract 

We propose a new method to compute the free energy or enthalpy of fluids or disordered 
solids by computer simulation . The main idea is to construct a reference system by freezing 
one representative configuration, and then carry out a thermodynamic integration. We present a 
strategy and an algorithm which allows to sample the thermodynamic integration path even in 
the case of liquids, despite the fact that the particles can diffuse freely through the system. The 
method is described in detail and illustrated with applications to hard sphere fluids and solids 
with mobile defects. 
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1. Introduction 

The free energy or free enthalpy is a central quantity in the statistical physics and thermody- 
namics of equilibrium systems, and there has traditionally been large interest in computing free 
energies in many areas of science [1|. Free energies need to be evaluated if one wishes to lo- 
cate first order phase boundaries accurately, to compare differently folded conformations of large 
molecules, or to characterize adsorption properties. In all of these cases, the quantities of interest 
are actually free energy differences, and a large number of elegant methods have been devised 
that allow to relate the free energies of specific systems with each other ||2]|3]|4][5]|6]|7][8][9]. 
An obvious alternative is to compute the absolute free energies of the two systems separately. 
Much less effort has been devoted to this type of approach. Computing the absolute free energies 
of arbitrary disordered or fluid systems such as, e.g., lipid membranes, solid with highly mobile 
defects, or nematic liquid crystals, remains one of the long-standing unsolved problems in com- 
puter simulations. If a fluid state has a direct connection with an ideal gas state in phase space, its 
free energy can be evaluated by standard thermodynamic integration techniques. In the systems 
mentioned above, however, the fluid phase of interest cannot be accessed from the gas phase in 
an obvious manner without crossing a first order phase boundary. 

Here we propose a method to calculate the free energies of arbitrary liquids and disordered 
solids [10 1, which will hopefully contribute to solve the problem. The main idea is to freeze 
an arbitrary configuration in the phase of interest - a 'representative' configuration, obtained, 
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e.g., from a typical simulation of an equilibrated system - and to use this frozen 'representative' 
configuration as the reference system for a thermodynamic integration. This idea seems bold 
and precautions have to be taken to make it work for fluids. We propose an algorithm that 
allows to determine the free energies for simple fluids and solids (as we shall demonstrate for 
hard spheres), which can presumably be generalized to complex fluids in a fairly straightforward 
manner. 

To put the method into context, we briefly sketch a few similar proposals in the past. Our 
idea is based on the "Einstein crystal" method originally developed by Frenkel and Ladd IfTTI . 
which allows one to calculate the free energy of a regular crystal structure by establishing a 
thermodynamic integration path to a reference system where particles are pinned to the lattice 
sites of the same crystal by harmonic springs. This method has been generalized in various ways. 
For example, the "Einstein molecule" |[T2l method alleviates the practical difficulties associated 
with shifts of the whole crystal during a simulation, and an "Einstein well" method was designed 
for "cluster crystals" where the lattice sites are occupied by several particles Ifl3ll . Speedy has 
generalized the Einstein method for amorphous solids by proposing to use reference systems 
with irregularly distributed sites [14], much in the spirit of our 'representative' reference system. 

The basic problem with applying the Einstein method to fluids is that "Einstein" particles are 
bound to their reference site - even in the limit of very weak springs - whereas fluid particles ex- 
plore the whole space. In infinite space, the transition from the "Einstein" system to the free fluid 
is thus singular. In finite space, the singularity disappears, but sampling the relevant quantities 
for the thermodynamic integration remains difficult (see below). Tyka et al lfl"5l have recently 
proposed to remedy this situation by using a complicated multi-particle spring potential, where 
particles are assigned to reference sites such that the total spring energy is minimal, under the 
constraint that there is still a one-to-one correspondence between particles and reference sites. 
Thus particles swap reference sites as they move along. As an application example, Tyka et al 
have calculated the free energies of liquid water and liquid argon. Their method has some sim- 
ilarities with the method we shall propose below, but our approach is conceptually simpler, and 
we do not need to solve a complex assignment problem for every single configuration. 

Our paper is organized as follows. In the next two sections, we present first the basic strategy 
of our method, and then the practical implementation. Extensions of the method, e.g., to constant 
pressure pressure simulations, are discussed in Section [3] In Section [4j we present first results 
regarding the behavior of the algorithm for different algorithmic parameters, and applications to 
hard sphere fluids. We conclude and summarize in section[5] 

2. Basic Strategy 

For clarity, we will constrain the discussion to monatomic liquids at constant volume at first. 
Extensions, e.g., to constant pressure simulations will be introduced further below (Section [3J. 
Furthermore, we only consider the configurational part of the free energy, since the kinetic part 
can be evaluated trivially [16]). The configurational energy is then characterized by a potential 
function or 'Hamiltonian' £/({r,)}) + H + Q, which depends on the coordinates {r, } of the particles 
i, and the configurational free energy shall be denoted Fq. 

The first step is to choose a 'representative' reference configuration, {r' ef }, which is typi- 
cally some arbitrary equilibrated configuration. The reference system is defined by the reference 
Hamiltonian 




(1) 
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where the function <S>(x) defines an attractive potential well with cD(x) < for x < 1 and O = 
elsewhere. The last term £/({r ref }) denotes the potential energy of hypothetical particles pinned 
to the reference sites and is thus constant. It is introduced for convenience. Note that there is 
formally a one-to-one correspondence between particles and reference sites, like in the Einstein 
crystal method. However, we want the particles to be indistinguishable (in the 'classical' sense), 
hence we allow them to swap identities (i.e., labels i and j) during the simulation. This indis- 
tinguishability together with the finite range of the well potential (in the Einstein crystal, it is 
infinite) helps to avoid the conceptual problems with the transition between the reference sys- 
tem H K f and the target system Ho mentioned earlier. Moreover, we shall show below (Section 



2.3 1 that the swap moves also speed up the equilibration times considerably. In the reference 



system, the particles do not interact directly with each other, hence its free energy F re f(s) can be 



calculated analytically (see Section 2.1 



To establish the connection between the reference system H Te f and the target system Ho, we 
introduce an intermediate model 

H'(A,s) = H inta .(A) + H Kf ( S ) (2) 

that reduces to the target system Ho at {A — 1,6 = 0) and to the reference system ff re f (e) at A = 
0. For example, one can choose H inter (A) = A AU with AU := t/({r,}) - U({rf})]. Alternatively, 
a more complicated function like HmtexW = -kgT ln(l - A + Ae~& ) can be more convenient 
when dealing with hard core interactions. The free energy of the intermediate system |2]) shall be 
denoted F'(A, s). Within the intermediate model, we can construct a thermodynamic integration 
path connecting the reference and target system. In practice, it is often simplest to evaluate 
separately the free energy difference associated with switching on the particle interactions (A = 
— > 1) at fixed s, and the free energy difference associated with switching off the well potentials 
in full presence of particle interactions (A - 1). Other integration paths are also conceivable and 
may be more effective in some cases. 

We will now discuss separately the different practical steps involved in this program. 

2.1. Free energy of the reference system 

We begin with calculating the free energy of the reference system. The partition function of 
a single particle in the volume V subject to a well potential <L>(|r - r ref [/r cuto ff) is given by 

Qi = V+ j dr [ exp(-/3e 0>) - l] =: V + V g<p(/te) 

with (3-1 /k B T, where Vq is the volume of a sphere with radius r cuto( f . This defines the function 



iM.Un =d f dn"- 1 [e-™ w - 1] (3) 
Jo 

characterizing the effect of a well potential <D(ji:) in d spatial dimensions on the single-particle 
partition function. Below, we shall mostly use a linear well potential, 

, . . _ f x - 1 for x < 1 
X) ~ \ otherwise, ' ' 



giving 

_ d! 



d k 

- a" 



k=0 kl 



(5) 
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In terms of the characterizing function g®(a), the configurational free energy of the reference 
system of N particles subject to the Hamiltonian H re f(s) finally reads 

+ t/({i? f }). (6) 

2.2. Switching on the particle interactions 

The next step consists in switching to the intermediate model H'(A, s) and evaluating the free 
energy difference associated with taking the interaction parameter A from zero to one. This could 
be done by a suitable thermodynamic integration. Here we will describe a different approach, 
which is also commonly used in the literature: We carry out a series of extended ensemble 
simulations where the parameter A can switch between two values Aq and A\ {e.g., in the sim- 
plest case Aq — and A\ = 1) by means of Monte Carlo trial moves, which are accepted or 
rejected according to a Metropolis criterion with the Hamiltonian H'(A,s) (Equation Q). The 
free energy difference of two systems with interactions Hi ntel (Ao) and H mta -(Ai) is then given by 
AF = F'(Ai,e) - F'(Aq,s) = -kgT ]n(P^JP^ ), where P^ is the fraction of configurations with 
interaction parameter /I,-. 

The question remains how to choose the steps A,. On the one hand, the free energy difference 
AF should not exceed several ksT, otherwise the system stays in the more favorable state forever. 
On the other hand, it should not be too small, otherwise it becomes difficult to control the relative 
statistical error. By a simple argument, we can estimate that the 'optimal' free energy difference 
should be about 2 kgT: 

Consider two systems with a total free energy difference AFo, and we have the task to evaluate 
the quantity q = exp(-AFo) by extended ensemble simulations of total length N. For simplicity, 
we assume that every Monte Carlo sweep results in an independent configuration. We compare 
two different simulation strategies: In the first, q is evaluated in one step within one single 
simulation run of length N, where the two systems are compared directly. In the second, q is 
evaluated in k steps, passing (k - 1) intermediate system with free energies equally distributed 
between the two systems of interest. This involves k simulation runs of length N/k. In the first 
setup (one-step simulation), the fraction of configurations N] in the state of higher free energy 
is given by NJN ~ {NJ(N - Ni)) = exp(-(3AF ) with the relative error ANJNi ~ VT/Wi" 
(assuming Poisson statistics). Thus the simulation yields the quantity q = exp(-fiAFo) with the 

relative error 

Aq/q = expOSAFo/2) • yfl/N. (7) 

In the second setup, (A:-step simulation), each simulation i (i = l,-,k) gives a quantity q t = 
exp(-/?AF 1 ) with the relative error Aqj/q, = exp(y8AF,/2) • y/k/N. Ideally, AFj = AF a /k for all i, 
and the total relative error of q — F]j <?; accumulates to 

Aq/q = exp(/3AF /2k) ■ k/ VA 7 . (8) 

Comparing ([7]) and (JHJ, one finds that a multi-step strategy is favorable as soon as the free energy 
difference AFo exceeds 41n(2) ~ 2.77 kgT, and the optimal step size is characterized by a free 
energy difference of f3AFo/k = 2k B T. 

2.3. Switching off the well potentials 

Finally, in the last step, we gradually turn the well potentials to zero. We evaluate the free 
energy difference AF2 between the intermediate model H'(A = l,s) and the target model H'(A = 



F Kf (s) = Nk B T 
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Figure 1: Sketch of moves in our Monte Carlo algorithm, (a) Simple particle displacements. (Could be replaced e.g., by 
short Molecular Dynamics runs.) (b) Smart particle swaps, (c) Smart particle relocations. See text for explanation 



1 , s — 0) = Hq by a thermodynamic integration, using 

Thus we need to sample the quantity (£, 0(|r, -r ref |/r cuto ff)) on an integration path between s' = s 
and s' = 0. At this point, it becomes clear why our potential wells must have a finite range r cato g : 
If the wells had infinite range, sampling (O) in fluids would be practically impossible in the limit 
s' = 0, where the mean-square displacement of the particles diverges. 

Setting a finite range for the reference potential, however, introduces a problem: The particles 
need to find their respective wells of attraction. We therefore introduce two Monte Carlo moves 
that help particles i explore their well i (Fig. [TJ. The first move (Fig. [T]b) swaps particles in a 
smart way and works as follows: 

• Pick a random particle i and find the set of particles {«, } that are within the attraction range 
of well i. 

• If particle i £ {«,}: pick a particle j from {«,}, and swap i and j with the probability 

min{l,£ e - A "'}. 

• Otherwise: pick a particle j from all particles 

- if j t {«,}: swap with probability min{l, ^e~ AH '}. 

- if j £ {«;}: swap with probability min{l, e~ AH }. 

Here AH' is the difference of the energies (according to the intermediate model) of the old and 
new configuration. This algorithm promotes particle swaps that bring particles close to their 
respective well and nevertheless satisfies detailed balance. 

The second move (Fig.[T]c) relocates panicles i with a bias towards the neighborhood of their 
well v. 

• Pick a random particle i (with position r,). 

• Choose a new position r' from a given (biased) distribution P,(r0 = exp(-W(|r' - rf |)). 

• Relocate the particle from r to rj with probability min{ 1 , P(r ; ) / P(rJ ) e~ AH ' } . 
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Figure 2: Illustration of the effect of the moves of Fig.[T]on the equilibration of a system of 80 hard disks (diameter D) 
at a density p = 0.8/D 2 , after switching on linear well potentials with strength e = SOk^T (r euto ff = 2D). Swap moves 
and relocation moves (one per bead) were attempted one per 100 MC sweeps. Left: Evolution of (<D) in simulations 
that include different moves as indicated. Right: Corresponding final configurations. Circles indicate particle positions, 
crosses give well positions. Particles and their respective wells are connected by straight lines. See text for more 
explanation. From Ref. 1101 . 

Obvious choices for W(r) which we have tested are W(r) = s(S>(r/r cuto s), or W(r) - const, for 
r < r cuto ff . As long as the wells are weak (low e'), the relocation moves have little effect. At high 
s', they help to overcome trapped situations where most particles are bound to a well, and a few 
cannot escape from a local cage. 

The effect of the different moves is shown in Fig. [2] (taken from Ref. |10|). It shows the 
evolution of the observable (<&,), averaged over all particles i, in a two dimensional system of hard 
disks (diameter D) at packing fraction r\ = 0.63, after s had been raised from zero to s = 50ksT. 
The well potential is linear (Eq. Q) with range r cutoi f = 2D. In a Monte Carlo simulation that 
includes only random particle displacements, the system is still far from equilibration after one 
million MC sweeps (a). Smart swap moves (one per 100 Monte Carlo sweeps) speed up the 
equilibration by orders of magnitude, but the system gets trapped in a configuration where one 
particle remains separated from its well (b). This problem is solved by including smart relocation 
moves (c). 
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3. Extensions of the method 

The algorithm described above can easily be generalized to molecular fluids: one must sim- 
ply swap molecules instead of particles (atoms). To generalize it to free enthalpy calculations 
(systems at constant pressure), we employ a reference system that is defined in terms of scalable 
coordinates, i.e., the reference sites are rescaled along with the particle coordinates if the vol- 
ume of the system changes. Furthermore, the volume V of the system is pinned to the reference 
volume by an additional term s(V - V re f) 2 in the reference Hamiltonian. We can then follow 
the strategy described above, with two modifications: First, the derivative of the free energy in 
Equation (|9|l has an additional contribution 

Second, the expression for the free enthalpy of the reference system changes. In practice, the 
fluctuations of the volume V become negligible already at moderate s, (e > IOIcbT in our case), 
such that the reference free energy can be written as 



N N \V re f/ \ V re f 

to a very good approximation. 



4. Results 

4.1. Transition to the reference state: Effect of well shape and well range 

Before presenting our first real application examples, we will discuss a few properties of the 
algorithm at the example of the two-dimensional ideal gas. Figure [3] illustrates the influence of 
the shape and range of the well potential. Three well shapes are compared, the (most popular) 
harmonic well, the linear well (Equation Q, and a "square root" well with the shape O(x) = 
yfx - \. The top panel in Figure [^demonstrates that of these three potentials, the harmonic well 
is most effective in trapping the particles for small and moderate e. All three potentials however 
bind almost 100 % of all particles for e > lOksT . At higher s, the linear and the square root 
potential localize the particles much better (lowest panel in Figure B), which makes it more easy 
to switch on interactions (step |2.2| > in the case of interacting particles. The middle panel in Figure 
[3] shows the averaged well energy per particle (0)/7V, i.e., the quantity that needs to be sampled 
for the thermodynamic interactions. It is remarkable that (<t>}/7V is still noticeably different from 
its limiting value ({<$>)/N = 1) even for e-values as high as s = 50IcbT. Despite the fact that 
it localizes particles most efficiently, the square well potential reaches the limiting value more 
slowly than the other potentials. In the following studies, we will mostly settle on using linear 
well potentials. 
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Figure 3: Effect of well shape Q>(x) on the behavior of different quantities as a function of well strength e in the two 
dimensional ideal gas. Top: Fraction / we u of particles inside the well. Middle: Free energy derivative per particle 
(®)/N. = -dF/ds ■ N . Bottom: Mean-square displacement of those particles that are within the well range. Dashed 
line corresponds to a quadratic well, thick solid line to a linear well, and dot-dashed line to a square root well with well 
range <x. The thin solid line shows the behavior of the same quantities for a linear well with range 2<t for comparison. 
The volume is V = 77.94cr 2 (Here tr is an arbitrary length scale). 
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4.2. Transition to the reference state: Effect of packing fraction 

Next we study the effect of particle interactions. To this end, we introduce hard disk inter- 
actions in the system of Figure [3] (using the linear well potential with range r cuto jf = 2<x). By 
adjusting the disk radius, we vary the packing fraction rj between r\ — and rj = 0.9, which is 
close to the maximum value ?7 max = 0.9069. At the two higher packing fractions, rj = 0.7 and 
77 = 0.9, the stable state is an ordered crystal structure. At rj — 0.7, we also show for comparison 
the curve obtained with a (metastable) disordered reference state. 

With increasing packing fraction, the particles settle more rapidly in their respective wells. 
This also happens faster for the ordered system than for the jammed system. To check for hys- 
teresis effects, the curves were sampled in both directions, with i.e., increasing and decreasing 
e. We did not encounter any signs of a first order transition. At the highest packing fraction 
77 = 0.9, the equilibration times were found to become very large for small e. The reason is that 
in the absence of wells, the center of gravity of the whole crystalline system diffuses around very 
slowly. This problem is already well known from the 'Einstein crystal' method and has led to the 
development of the 'Einstein molecule' method, where the reference system is defined in terms 
of relative coordinates rather than absolute coordinates |[T2l . The same idea can also be used 
here. 
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4.3. Application: Hard sphere fluids 

After these more technical discussions, we turn to presenting real applications. As a first ap- 
plication example, we have determined the absolute free energies of systems of 256 hard spheres 
with diameter D at the densities N/V = 0.25D 3 , 0.5D 3 , and 0.75Z)~ 3 . The results can be com- 
pared with those obtained by integration of the Carnahan-Starling equation-of-state[ 17 1, which 
describes the behavior of three dimensional hard sphere fluids very accurately. To obtain the 
free energies, we have used 50 values of e and 6 x 10 5 Monte Carlo sweeps for each data point 
for the two lower densities, and 200 values of e times one million Monte Carlo sweeps for the 
density 0.75D~ 3 . The results are given in Table [l] Within the error, all results agree well with 
the theoretical value from the Carnahan-Starling equation. 

For the density 0.5D -3 (corresponding to a liquid state), we compared three different situa- 
tions: (a) linear well <t> and liquid reference state, (b) linear well <t> and crystalline reference state, 
(b) harmonic well <t> and liquid reference state. All variations gave the same results. In case (b), 
we did not see a hysteresis on increasing/decreasing s - apparently, the trapping of particles in a 
crystalline array of wells is not associated with a phase transition at this density. We expect that 
this will be different closer to the liquid/solid transition. Nevertheless, we can conclude that our 
method is robust and may work if the reference configuration is not really representative for the 
target structure. Comparing the calculation for linear and harmonic wells, we find that the result 
also does not depend on the shape of the well. However, for more accurate potential, the linear 
potential seems to be more useful, because the particles get trapped more efficiently. 

Table 1: Results for the free energy of hard spheres (from Ref. 1101 1. F/N cs is the value according to the Carnahan- 
Starling equation-of-state [IV]. a) linear potential d>, liquid reference state, b) linear <I>, hep reference state, c) harmonic 
<1>, liquid reference state. 



N/V 


F/N 


(F/N) cs 


0.25 


0.620 ± 0.002 


0.625 


0.5 a) 


1.541 ±0.002 


1.544 


0.5"> 


1.540 + 0.002 


1.544 


0.5 C) 


1.549 ±0.002 


1.544 


0.75 


3.009 ± 0.005 


3.005 



4.4. Free energy of a vacancy 

Our last example shows an application of the method to a dense disordered system, where 
the dynamics is driven by cooperative processes. We have studied hard disks (two dimensions, 
diameter D) up to densities where the equilibrium state is solid, and enforced a vacancy defect 
by taking one particle out of an otherwise ordered configuration. To calculate the free energy 
of the vacancy, we must compare the enthalpy of this system with that of an ordered system at 
the same pressure. Thus simulations were carried out a constant pressure (see Section [3} in a 
rectangular simulation box of varying area and fixed side ratio 1 : V3/4. The vacancy is then 
stable, but highly mobile (Figure [5] top right). 
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Figure 5: Free enthalpy calculations of dense hard-disk systems (diameter D). Top left: Pressure vs. packing fraction 
as obtained from constant pressure simulations at e = 0. (a) expanded from an ordered solid phase, (b) expanded from 
a solid phase with one vacancy, (c) compressed from the fluid phase. The solid line shows the theoretical estimate [ 1 8 1 
P = p/(l - p ■ 7r/4) 2 with p = (N + 1)/(V). Top right: Configuration snapshots with one vacancy at the beginning 
(crosses) and the end (circles) of a Monte Carlo run at e = 0. The arrows mark the position of the defect in the two cases, 
illustrating its high mobility. Bottom: Derivative -dF/de of the free energy along the path of thermodynamic integration 
for the three dense systems (a),(b), and (c) at pressure P = WksT/D 2 , and for a fluid system at P = 6kgT /D 2 . The inset 
shows the same data in a semi-logarithmic scale. The circle highlights the area where the volume fluctuations contribute 
noticeably to dF/de. 
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We compare four different structures (Figure |5J: An ordered solid (a), an ordered solid with 
a vacancy (b), and for comparison also a metastable disordered jammed state (c), which was 
obtained by compressing the system from the fluid state, and a fluid state (d) at lower pressure. 
Figure|5](top left) shows the corresponding pressure-density curves of the target system (s = 0). 
Free enthalpy calculations were carried out at P = l0k B T/D 2 for the three cases (a-c), and at 
P = 6k B T/D 2 for the case (d). The resulting free enthalpies can be related to the chemical 
potential fj. by virtue of the thermodynamic relation G - fj.N. Figure [5] shows the corresponding 
evolution of the integrands (dH'/ds) of the thermodynamic integration (Equation (10 1 up to 
s = 20. (in total, s was varied between and 50k B T). The volume fluctuations contribute only 
for very small s. The curves for the jam, solid, and defect state are slightly different, reflecting 
the differences in the free enthalpies of the different states. 

At P — 6k B T/D 2 , the free enthalpy calculation yields the chemical potential yt = 8.997 ± 
Q.QQ2k B T, which is in good agreement with the theoretical estimate of Helfand et al |[T8l n = 
9.041k B T. At P = l0k B T/D 2 , we find // 80 ud = 13.617 + 0.002k B T in the solid state, and <i iam = 
13.675 ± 0.002k B T in the jammed state, which establishes that the solid state is indeed the stable 
phase. For the system with one defect, we obtained the total enthalpy Gdefect = 1361.7 + Q.2k B T . 
This result can be used to estimate the core free energy of the vacancy fi c - Gdefect - HsoiidN + 
k B Tln(N) = 7.1 ± Q3k B T, which corresponds to a relative vacancy frequency of roughly 10~ 3 . 
(For comparison, the frequency of vacancies at liquid/solid coexistence in three dimensions is 
roughly 10~ 4 , according to Pronk et q fl[T9"l ). Since finite size effects stabilize the defect, /i c is 
probably largely overestimated, hence the value given above has to be interpreted as an upper 
bound. More detailed studies shall be carried out in the future. Here, the example mainly serves 
to illustrate the use of our approach in situations where free energies are difficult to access with 
other methods. 



5. Conclusions 

In summary, we propose a new method to compute absolute free energies for fluids, liquids, 
and disordered structures. The method generalizes the Einstein crystal method of Frenkel and 
Ladd ifTTI . but can be applied much more widely. We anticipate that the method will be useful 
for studies of complex fluids and biological systems, e.g., proteins in solution. This remains to 
be tested. The method has recently been published IfTOll . and a few application examples have 
been given. We hope that the detailed discussion in the present proceedings paper will motivate 
researchers to try and apply the method to many other systems and further explore its potential. 
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